birds <- read.csv('ecol 563/birds.csv') birds.short <- birds[!is.na(birds$S),] # create variables for WinBUGS y <- birds.short$S year2 <- as.numeric(birds.short$year==2006) year3 <- as.numeric(birds.short$year==2007) n <- length(y) # number of patches J <- length(unique(birds.short$patch)) # patch identifier renumbered patch <- as.numeric(factor(birds.short$patch)) setwd("C:/Users/jmweiss/Documents/ecol 563/models") ##### BUGS model: model2.txt #### model{ for(i in 1:n) { y[i]~dpois(mu.hat[i]) log.mu[i]<-a[patch[i]] + b1*year2[i] + b2*year3[i] mu.hat[i] <- exp(log.mu[i]) } for(j in 1:J){ a[j]~dnorm(0,.000001) } b1~dnorm(0,.000001) b2~dnorm(0,.000001) } ############### bird.data <- list("y", "year2","year3","n", "J", "patch") bird.inits <- function() {list(a=rnorm(J), b1=rnorm(1), b2=rnorm(1))} bird.parms <- c("a", "b1", "b2") # WinBUGS run library(arm) mod1.bugs <- bugs(bird.data, bird.inits, bird.parms, "model2.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=10000, debug=T) # JAGS run library(R2jags) mod1.jags <- jags(bird.data, bird.inits, bird.parms, "model2.txt", n.chains=3, n.iter=1000) names(mod1.bugs) # WinBUGS objects returned by JAGS names(mod1.jags$BUGSoutput) # examine traces of Markov chains mod1.mcmc <- as.mcmc.list(mod1.bugs) xyplot(mod1.mcmc, layout=c(6,6)) mod1jags.mcmc <- as.mcmc(mod1.jags) xyplot(mod1jags.mcmc, layout=c(5,5)) # random intercepts model--first parameterization ########### BUGS model: model3.txt ########### model{ for(i in 1:n) { y[i]~dpois(mu.hat[i]) log.mu[i]<-a[patch[i]] + b1*year2[i] + b2*year3[i] mu.hat[i] <- exp(log.mu[i]) } for(j in 1:J){ a[j]~dnorm(mu.a,tau.a) } b1~dnorm(0,.000001) b2~dnorm(0,.000001) mu.a~dnorm(0,.000001) tau.a <- pow(sigma.a,-2) sigma.a~dunif(0,10000) } ##################### bird.inits <- function() {list(a=rnorm(J), b1=rnorm(1), b2=rnorm(1), mu.a=rnorm(1), sigma.a=runif(1))} bird.parms <- c("a", "b1", "b2", "mu.a", "sigma.a") mod3.bugs <- bugs(bird.data, bird.inits, bird.parms, "model3.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=50) mod3.bugs <- bugs(bird.data, bird.inits, bird.parms, "model3.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=1000, debug=T) max(mod3.bugs$summary[,"Rhat"]) min(mod3.bugs$summary[,"n.eff"]) mod3.bugs$summary[101:106,] # random intercepts model--second parameterization ##### BUGS model: model3a.txt ###### model{ for(i in 1:n) { y[i]~dpois(mu.hat[i]) log.mu[i]<-a[patch[i]] + b1*year2[i] + b2*year3[i] mu.hat[i] <- exp(log.mu[i]) } # frequentist parallel--random intercept effects for(j in 1:J){ u0[j]~dnorm(0,tau.a) a[j] <- mu.a + u0[j] } b1~dnorm(0,.000001) b2~dnorm(0,.000001) mu.a~dnorm(0,.000001) tau.a <- pow(sigma.a,-2) sigma.a~dunif(0,1000) } ############### bird.inits <- function() {list(u0=rnorm(J), b1=rnorm(1), b2=rnorm(1), mu.a=rnorm(1), sigma.a=runif(1))} bird.parms <- c("u0", "b1", "b2", "mu.a", "sigma.a") mod3a.bugs <- bugs(bird.data, bird.inits, bird.parms, "model3a.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=1000, debug=T) # frequentist version of the model library(lme4) mod.lmer <- lmer(S~factor(year)+(1|patch), data=birds, family=poisson) fixef(mod.lmer) # frequentist patch random effects round(ranef(mod.lmer)$patch[,1],3) # compare with Bayesian estimates round(mod3a.bugs$summary[,"mean"],3) # frequentist random intercepts coef(mod.lmer ### add patch-level predictors: area and landscape ### unique(birds$landscape) mod2.lmer <- lmer(S~factor(year)+log(area)+landscape+(1|patch),family=poisson,data=birds) # obtain one value per patch L.area <- tapply(log(birds.short$area),birds.short$patch, function(x) x[1]) # L.area is an array, we need it to be a vector class(L.area) L.area <- as.numeric(tapply(log(birds.short$area), birds.short$patch, function(x) x[1])) class(L.area) # there is a missing value due to a phantom factor level L.area # redeclare patch as a factor L.area <- as.numeric(tapply(log(birds.short$area), factor(birds.short$patch), function(x) x[1])) # create landscape variable Lscape <- as.numeric(tapply(birds.short$landscape, factor(birds.short$patch), function(x) x[1])) # dummy variables for landscape Lscape2 <- as.numeric(Lscape==2) Lscape3 <- as.numeric(Lscape==3) Lscape4 <- as.numeric(Lscape==4) #### BUGS program: model4.txt ########### model{ for(i in 1:n) { y[i]~dpois(mu.hat[i]) log.mu[i]<-a[patch[i]] + b1*year2[i] + b2*year3[i] mu.hat[i] <- exp(log.mu[i]) } for(j in 1:J){ a[j]~dnorm(a.hat[j],tau.a) a.hat[j] <- mu.a + g1*L.area[j] + g2*Lscape2[j] + g3*Lscape3[j] + g4*Lscape4[j] } g1~dnorm(0,.000001) g2~dnorm(0,.000001) g3~dnorm(0,.000001) g4~dnorm(0,.000001) b1~dnorm(0,.000001) b2~dnorm(0,.000001) mu.a~dnorm(0,.000001) tau.a <- pow(sigma.a,-2) sigma.a~dunif(0,10000) } #################### # change data, inits, and parm objects bird.data <- list("y", "year2", "year3","n", "J", "patch", "L.area", "Lscape2", "Lscape3", "Lscape4") bird.inits <- function() {list(a=rnorm(J), b1=rnorm(1), b2=rnorm(1), mu.a=rnorm(1), sigma.a=runif(1), g1=rnorm(1), g2=rnorm(1), g3=rnorm(1), g4=rnorm(1))} bird.parms <- c("a", "b1", "b2", "mu.a", "sigma.a", "g1", "g2","g3","g4") mod3a.bugs <- bugs(bird.data, bird.inits, bird.parms, "model4.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=100, debug=T) mod3a.bugs <- bugs(bird.data, bird.inits, bird.parms, "model4.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=1000, debug=T) max(mod3a.bugs$summary[,"Rhat"]) min(mod3a.bugs$summary[,"n.eff"]) mod3a.bugs$summary[100:110,] #compare to frequentist estimates fixef(mod2.lmer) mod3a.bugs$summary[102:110,] # calculate mean richness by landscape for 2005 and patches of area = 1 ##### alter BUGS program to calculate means: model4.txt ##### model{ for(i in 1:n) { y[i]~dpois(mu.hat[i]) log.mu[i]<-a[patch[i]] + b1*year2[i] + b2*year3[i] mu.hat[i] <- exp(log.mu[i]) } for(j in 1:J){ a[j]~dnorm(a.hat[j],tau.a) a.hat[j] <- mu.a + g1*L.area[j] + g2*Lscape2[j] + g3*Lscape3[j] + g4*Lscape4[j] } g1~dnorm(0,.000001) g2~dnorm(0,.000001) g3~dnorm(0,.000001) g4~dnorm(0,.000001) b1~dnorm(0,.000001) b2~dnorm(0,.000001) mu.a~dnorm(0,.000001) tau.a <- pow(sigma.a,-2) sigma.a~dunif(0,10000) # landscape means in 2005 when area=1 mean1 <- exp(mu.a) mean2 <- exp(mu.a+g2) mean3 <- exp(mu.a+g3) mean4 <- exp(mu.a+g4) } #################### bird.parms2 <- c("a", "b1", "b2", "mu.a", "sigma.a", "g1", "g2","g3","g4", "mean1","mean2","mean3","mean4") # rerun model using last values of previous run mod4.bugs <- bugs(bird.data,mod3a.bugs$last.values , bird.parms2, "model4.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=1000, debug=T) # percentile credible intervals apply(mod4.bugs$sims.matrix[,c("mean1","mean2","mean3","mean4")],2,function(x) quantile(x,c(.025,.975))) # HPD credible intervals HPDinterval(as.mcmc(mod4.bugs$sims.matrix[,c("mean1","mean2","mean3","mean4")])) # we can also just add the appropriate columns of sims.matrix # and use quantile function on the result quantile(exp(mod4.bugs$sims.matrix[,"mu.a"]+mod4.bugs$sims.matrix[,"g2"]),c(.025,.975))